Single-particle and collective slow dynamics of colloids in porous confinement 
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Using molecular dynamics simulations we study the slow dynamics of a hard sphere fluid confined 
in a disordered porous matrix. The presence of both discontinuous and continuous glass transitions 
as well as the complex interplay between single-particle and collective dynamics are well captured 
by a recent extension of mode-coupling theory for fluids in porous media. The degree of universality 
of the mode-coupling theory predictions for related models of colloids is studied by introducing 
size-disparity between fluid and matrix particles, as well as softness in the interactions. 
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Investigations of the dynamic properties of fluids con- 
fined in a disordered porous matrix have recently at- 
tracted considerable attention in the scientific commu- 
nity both from a fundamental and from a technological 
point of view. Of particular interest is the glass forma- 
tion process where the single-particle and the collective 
dynamics of the fluid are exposed to the complex inter- 
play between confinement and connectivity of the pores. 
While a considerable amount of experimental results has 
been compiled in this field [l], 0] , theoretical investiga- 
tions are still rare: remarkable computer simulations on 
a few selected state points and systems have been per- 
formed 0, Q , but systematic studies of this phenomenon 
are difficult due to the large computational cost caused 
by the requirement of long simulations and by the ther- 
modynamic averaging prescriptions for these systems. 

A breakthrough in the theoretical investigations was 
achieved with the formalism put forward by Krakovi- 
ack [H @, 0] who successfully combined two concepts 
to study the dynamics of fluids confined in disordered 
porous matrices. On one hand, his framework is based on 
the replica Ornstein-Zernike (ROZ) formalism where 
the system at hand is viewed as a quenched-annealed 
(QA) mixture of mobile fluid particles and immobile ma- 
trix particles. The other ingredient is mode-coupling the- 
ory (MCT) Q: based on the static correlations of the 
system, MCT predicts the time dependence of the single- 
particle ("self") and collective density correlators. MCT 
has been able to identify new and unexpected features 

in the dynamics of colloidal systems [13, IHL E3, E3] and 

provides a convincing account of the early stages of the 
structural arrest in molecular glass-forming liquids 14 1. 

The results for the kinetic diagram of a simple hard 
sphere (HS) fluid in a HS matrix, evaluated using the 
MCT framework with ROZ structure factors as an input, 
contains indeed a wealth of intriguing features @, 0] : (i) 
two types of glass transitions are encountered, a discon- 
tinuous type B transition at low matrix packing fractions 
(</>,„), and a continuous type A transition at large <j> m [231 ]: 
the two transition lines meet at a degenerate high order 
singularity point (A3). Further, the theory predicts (ii) 
a re-entrant glass transition for large <f> m and (iii) a con- 



tinuous diffusion-localization transition that occurs only 
in the self dynamics; the latter transition is driven by 
the localization of particles in disconnected void domains 
formed by the matrix 0], as in the Lorentz model [l5j |. 

This Letter aims in two directions: on one side we 
provide for the first time a parallel study to Krakovi- 
ack's main theoretical predictions by means of computer 
simulations for the same HS system. The simulations 
provide, in addition, a "realistic" picture of the above 
dynamic effects as opposed to the "idealized" view that 
emerges from MCT. On the other hand, with the help 
of the ROZ+MCT framework we explore the parameter 
space of two further types of potentials to find out to 
which degree the predicted dynamic features are generic. 
To this end we mimic typical interactions encountered in 
colloidal systems by introducing a size disparity between 
matrix and fluid particles (establishing thereby a link to 
the size- asymmetric mixtures investigated in (l6l.[l7. 18 1) 
and softness in the interparticle interaction. 

By means of event- driven molecular dynamics (MD) we 
study a QA system of hard spheres in three dimensions 
under periodic boundary conditions. As in the monodis- 
perse model studied in [fj , both matrix and fluid particles 
have the same diameter a rn — <jf — a and mass m. The 
total number of particles was N = Nf + N m = 1000. In 
the following, we use a and a/ 'fcsT '/ 'ma 2 as units of length 
and time, respectively (l9| . The porous matrix was gen- 
erated by quenching an equilibrium fluid configuration 
at a packing fraction <f> m — ir j^>p m a^ n . Subsequently, the 
fluid component with packing fraction cf>f — n/Gpfaj: was 
inserted into the remaining volume using an optimized 
inflation algorithm. For each state point, dynamic prop- 
erties were averaged over ten independent realizations of 
the system. The fluid component was considered equi- 
librated if the total mean squared displacement (MSD) 
exceeded 10 2 within a simulation time of 3 x 10 4 . State 
points for which at least half of the system realizations 
did not meet this criterion were considered glassy. No 
signs of crystallization were observed. The MCT equa- 
tions for QA systems were solved using the numerical 
procedure outlined in @, 0] . For the HS system, the con- 
nected and blocked structure factors, used as an input 
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FIG. 1: Kinetic diagram of the OA HS system: MCT tran- 
sition lines (thin lines; from Ref. Q) and dynamic arrest line 
from simulations (thick line; see text for definition). The A3 
singularity predicted by MCT is indicated by a filled trian- 
gle. Empty and filled circles indicate equilibrated and glassy 
simulated state points, respectively. 



for the MCT equations, were determined by solving the 
ROZ equations [§| with the Percus-Yevick (PY) closure 
relation. The so-obtained critical packing fractions and 
exponents typically agree within few percent with those 
reported in [fl 0] . 

In Fig. [TJ we show in the ((f) m , <fif) -plane the glass and 
the diffusion-localization transition lines obtained from 
MCT 0, 0] together with the dynamic arrest line de- 
termined from simulations. The arrest line interpolates 
through points with adjacent equilibrated and glassy 
state points. The shape of the arrest line resembles qual- 
itatively the ideal MCT glass line, but it starts bending 
downwards more rapidly for <p m > 0.2. Performing sim- 
ulations on a fine mesh in <pf and <f> m we found that the 
arrest line is not re-entrant at small <f>t. We will show 
below that this apparent discrepancy can be rationalized 
in terms of a crossover from glass to diffusion-localization 
transition in the self dynamics. 

To investigate in more detail the different transition 
scenarios predicted by MCT, we focus on selected paths 
(indicated by arrows and labelled in Fig. [J) across the 
(</>,„, 0/)-plane. We study both the self part F s (k,t) — 
1/AT(YJ . exp{ik- [r 3 -(t) — rj-(O)]}) and the connected part 
F c (k,t) = (W*MPk(°MPk(0)/o£(0)) of the intermedi- 
ate scattering function [6(, where 5pk(t) = Pk{t) — (pk), 
Pk(t) — J2j exp[ik • rj(t)], and (...) indicates a thermal 
average. F c (k,t) is the appropriate correlator to charac- 
terize the slow collective dynamics of fluids in porous me- 
dia, since the conventional coherent intermediate scatter- 
ing function, F(k,t), fails to decay to zero at long times 
due to matrix- induced average density fluctuations Q . 
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FIG. 2: Upper panels: intermediate scattering functions at 
4> m = 0.05 for different values of 4>f. (a) connected part 
F c (k = 7.0, t) and (b) self part F s (k = 7.0, t). Solid lines 
are fits to the second-order expansion of the /3-correlator (see 
text). Error bars for selected times represent one standard de- 
viation on the average over matrix realizations. Lower panels: 
non-ergodicity parameter f(k) (symbols) for (c) connected 
and (d) self correlators at cf> m — 0.05, <j>f — 0.50. The cor- 
responding non-ergodicity parameters predicted by MCT at 
4>f — 0.446 are included as solid lines. 



Figure [2] shows the evolution of the intermediate scat- 
tering functions upon increasing (j)f at fixed <p m = 0.05 
(path I). The wave-vector considered is k = 7.0, close to 
the main peak of the fluid-fluid static structure factor. 
Both the self and the collective dynamics slow down sig- 
nificantly as 4>f increases towards 0.51, while a plateau at 
intermediate t develops both in F s (k,t) and in F c (k 7 t). 
These features are similar to those observed in bulk glass- 
forming liquids and are compatible with a "smeared" 
type B transition. The difference between the MCT criti- 
cal packing fraction, <^ = 0.439645, and the glass transi- 
tion observed in our simulations amounts to w 15%. Sim- 
ilar discrepancies have been reported in investigations on 
bulk hard spheres and suspensions of hard colloids 



14] 



To test the predictions of MCT at a more quanti- 
tative level, we describe the dynamics within the /3- 
relaxation regime (close to the plateau) using the second- 
order asymptotic expansion of the MCT /3-correlator, 
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FIG. 3: Intermediate scattering functions at <f> m = 0.20 for 
different values of cf>f-. (a) connected part F c (k — 7.0, i) and 
(b) self part F s (k — 7.0, t). Error bars as in Fig. [2] Inset 
of (b): enlarged view of the plateau of F s (k = 7.0, t) at long 
times. 



f{k) + Hi(k)t b + H 2 (k)t 2b 9]. We perform a simulta- 
neous fit to all integer wave-vectors 1 < k < 20 at the 
largest equilibrated packing fraction, <f)f = 0.50. The re- 
sult for k = 7.0 is included in Fig. [21 the exponent b is 
0.33 ±0.05 for F a (k,t) and 0.36 ± 0.05 for F c (k,t). The 
range over which the fit can be considered reliable cov- 
ers « 1.5 decades in time. The estimated uncertainty on 
b accounts for the sensitivity of the fit to the choice of 
the time range. The fitted values are somewhat smaller 
than b — 0.51 (independent of correlators) obtained from 
numerical solution of the MCT equations, which may be 
due to the limited time window available for fitting the 
simulation data. Further, in the bottom panels of Fig. [2] 
we show the non-ergodicity parameter f(k) for both self 
and connected correlators. In both cases, f(k) compares 
well with the MCT predictions at </>/ = 0.446, i.e., only 
slightly above (fy'j. Comparisons of simila r q uality have 
been reported for bulk glass-formers [2(J HI]]. We con- 
clude that MCT is able to describe at a semi-quantitative 
level the slow dynamics of the QA system in this portion 
of the kinetic diagram. 

A more complex scenario appears at large <p m , where 
the stronger influence of the matrix structure changes the 
nature of the transition. In Fig. Owe show the intermedi- 
ate scattering functions for k = 7.0 at constant cf> m = 0.2 
(path II in Fig. [TJ) . Along this path, the relaxation pat- 
terns for self and collective correlators differ markedly as 
4>f is varied. The relaxation of F c (k,t) becomes slower 
as <j)f increases. Contrary to the case of dilute matrices, 
however, there is no finite-height plateau at intermedi- 
ate times. This indicates the approach to a continuous 
type A glass transition, consistent with MCT predictions 
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FIG. 4: MSD at constant <j> } = 0.1 for <j> m = 0.05, 0.10, 0.15, 
0.20, 0.225, 0.25, and 0.275 (from left to right). Solid lines are 
fits to 5r 2 (t) ~ t" in the sub-diffusive regime. Error bars as in 



Fig. [21 Inset: z as a function of 
as in the main plot). 



for selected 4> m (symbols 



along path II. On the other hand, the decay of F s (k,t) 
occurs in two steps, with a first inflection (1 < t < 10) 
and a subsequent stretched decay at long times to a finite 
plateau. Upon increasing cj>f the plateau rises continu- 
ously starting from values close to 0, as one expects in a 
type A transition. These features of F s (k,t) are consis- 
tent with those predicted by MCT across the continuous 
diffusion- localization transition (dotted line in Fig.[TJ). 

The cascade of decays in F s (k, t) hence arises from the 
interplay between two dynamic effects: a weak, collec- 
tive mechanism of caging by fluid particles and an ef- 
fective trapping in the voids left by the matrix particles. 
The superposition of glass and diffusion-localization tran- 
sitions leads to an effective "decoupling" between self 
and collective dynamics. In fact, at the largest equili- 
brated 4>f along path II ((f) f = 0.22) the ratio T s (k) / r c (k) 
of the relaxation times, defined via the decay to 0.1 
of the corresponding correlators, is already larger than 
10 2 for k — 7.0. As a result, the dynamic arrest line 
(see Fig. [TJ) bends downwards rapidly in the large-<^ m , 
small-0/ part of the diagram, following the trend of the 
diffusion-localization line. Further work is required to 
assess whether the predicted re-entrance of the glass line 
exists or if it is hindered by the progressive arrest of the 
self dynamics. 

According to MCT, the stretched decay of self corre- 
lators for large t close to the diffusion-localization line 
is associated with a sub-diffusive behavior in the MSD 
5r 2 (t) 0. In Fig. [4] we show the evolution of the MSD 
at constant <j)f — 0.10 (path III in Fig. [TJ). Note that 
upon increasing <f> m along this path, the slowing down of 
the self dynamics is dictated by the diffusion-localization 
transition. For <j) m < 0.15, the ballistic regime Sr 2 (t) ~ 
i 2 is followed by normal diffusion Sr 2 (t) ~ t. However, 
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FIG. 5: MCT kinetic diagram for (a) QA hard spheres with 
size ratio A = 0.5, 1.0, 2.0 (from left to right) and (b) QA soft 
spheres with A = 1.0 and repulsive exponent n — 6, 7, 10 (from 
top to bottom). Lines as in Fig[T] Filled triangles indicate 
points where type A and B transitions meet. 

for <f> m > 0.15, the system becomes sub-diffusive, i.e., 
5r 2 (t) ~ t z with z < 1, over an intermediate time win- 
dow that becomes broader as <p m increases. As in the 
Lorentz model [HI], z decreases as <f> m increases. Up 
to a (j)* n f=a 0.25, normal diffusion is eventually recov- 
ered at longer times, while for larger <f> m the MSD tends 
to saturate. Close to 4>* m: the sub-diffusion exponent is 
z « 0.51, in striking agreement with the MCT prediction 
along the diffusion-localization line [z = 0.5 [7]). Re- 
markably, we also find that z is nearly independent of <j)f 
at fixed (f> m (see inset of Fig. 01 . Hence, even at finite </)/, 
the diffusion-localization transition is intimately related 
to the geometrical properties of the matrix. 

To assess the degree of universality of the scenario 
predicted by the ROZ+MCT formalism, we solved the 
MCT equations for related QA systems based on purely 
repulsive interactions. Specifically, we considered addi- 
tive hard spheres with size ratio A = a m /af ^ 1 and 
soft spheres interacting with inverse power potentials 
u(r) = e(a/r) n where 6 < n < 12. In the latter case, we 
used A = 1 and the structure factors were calculated at 
a temperature T/e = 0.2 by solving the ROZ equations 
in combination with the hypernetted-chain approxima- 
tion [3] • fn Fig. [5] we show results for selected values 
of A and n. The topology of the MCT kinetic diagrams 
of the studied systems is very similar to the one shown 
in Fig. [TJ but subtle differences are observed as system 
parameters vary. As A decreases in the HS systems, the 
glass line shifts to smaller <f> m and the crossing point be- 
tween type B and A transitions moves towards the turn- 
ing point of the glass line. On the other hand, in the 
soft sphere systems the glass line shifts to larger 0/ as 
n decreases (as expected), and the diffusion- localization 
line becomes clearly re-entrant. We plan to investigate 
these features in more detail in future work. 

In conclusion, we revealed the existence of a subtle 
interplay between self and collective dynamics in a QA 
model of fluid in porous confinement. The MCT frame- 
work for QA systems 0, @] provides idealized but con- 



sistent guidelines to explain the dynamic features appar- 
ent for both dilute and dense matrices. This includes 
the superposition of glass transitions, driven by a collec- 
tive caging mechanism, and diffusion-localization transi- 
tions, associated with the ceasing percolation of voids in 
the matrix structure. The predicted features should be 
generic for a broad class of model colloidal fluids adsorbed 
in porous media, and should be observed experimentally 
in colloidal suspensions confined in porous matrix config- 
urations, quenched for instance by optical tweezers. 
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dation (FWF) under Proj. No. P19890-N16. 

Note — While finalizing this paper, we became aware 
of a recent paper by Kim et al. [22J , which focuses on a 
related model. Their results for the HS QA system are 
consistent with the analysis presented here. 
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